#####Load Libraries and data#####
library(pacman)
p_load(tm, ggplot2, SnowballC, wordcloud, topicmodels, stm, plyr, dplyr, 
       tidytext, quanteda,readtext, geometry, Rtsne, rsvd, erer, foreign, 
       ISLR, tidyr, magrittr, lubridate, caret, tictoc, scales, stats, purrr, 
       furrr, tictoc, mltools, Matrix, stminsights)

#####START HERE to Reload Data#####
#save.image("welfarePaper.RData")
load("welfarePaperFinal.RData")

#####Estimate Effects#####
set.seed(11042020)
prep <- estimateEffect(formula = ~ s(days) +partyControl+Party+age+
                         gender+hispanic+black+numberTerms+medianIncome+
                         prcntBA+prcntHS+prcntBlack+prcntHisp+
                         prcntWhite+prcntAsian+prcntOld+prcntUnemp+
                         recentArrivalPrcnt+totalEmploymentPop+region,
                       stmobj = stm74, 
                       metadata = meta)
beepr::beep(sound =2)
prep$data$days <- as.numeric(prep$data$days)

##bonferroni corrected p values for 74 topics##
summary.estimateEffect <- function(object, topics=NULL, nsim=500, ...) {
  if(is.null(topics)) topics <- object$topics
  if(any(!(topics %in% object$topics))) {
    stop("Some topics specified with the topics argument are not available in this estimateEffect object.")
  }
  tables <- vector(mode="list", length=length(topics))
  for(i in 1:length(topics)) {
    topic <- topics[i]
    sims <- lapply(object$parameters[[which(object$topics==topic)]], function(x) rmvnorm(nsim, x$est, x$vcov))
    sims <- do.call(rbind,sims)
    est<- colMeans(sims)
    se <- sqrt(apply(sims,2, stats::var))
    tval <- est/se
    rdf <- nrow(object$data) - length(est)
    p <- 2 * stats::pt(abs(tval), rdf, lower.tail = FALSE)
    p.adj <- ifelse(test = round(p*74,3) > 1, yes = 1,
                    no = round(p*74,3))
    
    coefficients <- cbind(est, se, tval, p,p.adj)
    rownames(coefficients) <- attr(object$parameters[[1]][[1]]$est, "names") 
    colnames(coefficients) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)","P.Adj")
    tables[[i]] <- coefficients
  }
  out <- list(call=object$call, topics=topics, tables=tables)
  class(out) <- "summary.estimateEffect"
  return(out)
}

#####Extraction Function for Plotting in GGPLOT#####
extract.estimateEffect <- function (x, covariate, model = NULL, topics = x$topics, method = "pointestimate", 
                                    cov.value1 = NULL, cov.value2 = NULL, moderator = NULL, moderator.value = NULL, 
                                    npoints = 100, nsims = 100, ci.level = 0.95, custom.labels = NULL, 
                                    labeltype = "numbers", n = 7, frexw = 0.5) 
  {
  cthis <- stm:::produce_cmatrix(prep = x, covariate = covariate, 
                                 method = method, cov.value1 = cov.value1, cov.value2 = cov.value2, 
                                 moderator = moderator, moderator.value = moderator.value)
  simbetas <- stm:::simBetas(parameters = x$parameters, nsims = nsims)
  uvals <- cthis$cdata[[covariate]]
  offset <- (1 - ci.level)/2
  labels <- stm:::createLabels(labeltype = labeltype, covariate = covariate, 
                               method = method, cdata = cthis$cdata, cov.value1 = cov.value1, 
                               cov.value2 = cov.value2, model = model, n = n, topics = x$topics, 
                               custom.labels = custom.labels, frexw = frexw)
  out <- lapply(topics, function(i) {
    sims <- cthis$cmatrix %*% t(simbetas[[which(x$topics == 
                                                  i)]])
    if (method == "difference") {
      diff <- sims[1, ] - sims[2, ]
      out_inner <- data.frame(method = method, topic = i, 
                              covariate = covariate, covariate.value = paste0(cov.value1, 
                                                                              "-", cov.value2), estimate = mean(diff), std.error = sd(diff), 
                              ci.level = ci.level, ci.lower = quantile(diff, 
                                                                       offset), ci.upper = quantile(diff, 1 - offset), 
                              label = labels[which(x$topics == i)])
    }
    else {
      #browser()
      out_inner <- data.frame(method = method, topic = i, 
                              covariate = covariate, covariate.value = uvals, 
                              estimate = rowMeans(sims), std.error = sqrt(apply(sims, 
                                                                                1, stats::var)), 
                              ci.level = ci.level, ci.lower = apply(sims,
                                                                    1, quantile, probs = offset), ci.upper = apply(sims, 
                                                                                                                   1, quantile, probs = (1 - offset)), label = labels[which(x$topics == 
                                                                                                                                                                              i)])
      out_inner["tval"] <- out_inner["estimate"]/out_inner["std.error"]
      rdf <- nrow(x$data) - length(out_inner["estimate"])
      out_inner["p"] <- (2 * stats::pt(abs(out_inner[["tval"]]), rdf, lower.tail = FALSE))
    }
    if (!is.null(moderator)) {
      out_inner$moderator <- moderator
      out_inner$moderator.value <- moderator.value
    }
    rownames(out_inner) <- NULL
    return(out_inner)
  })
  out <- do.call("rbind", out)
  return(out)
}


#####LABEL SUBSTANTIVE TOPICS#####
substantive_top <- c("Water Policy","Medicare","Equal Rights","Govt Spending","Conservation",
                     "Native Policies","Energy Policy","Agricultural Policy","War in Middle East",
                     "Veterans' Policy","Election Rules","Transportation","Healthcare","Homeland Security",
                     "Peace Negotiations","Appropriations","Legal Provisions","Self-Determination",
                     "Earmarks","Land Management","Drug Policy","Small Business","Trade Policy",
                     "Criminal Justice","Labor Policy","Abortion","Individual Taxation","Financial Policy",
                     "Human Rights Policy","Judicial Policy","Impeachment","Telecommunications",
                     "Technology Funding","Budget Deficit","Child Welfare","Medical Funding",
                     "Nuclear Policy")
substantive_num <- c(01,03,04,09,11,13,15,16,17,18,19,23,24,26,27,28,29,34,39,42,44,46,48,
                     49,53,54,55,56,58,59,61,62,63,65,68,71,73)

topics_and_nums <- data.frame(topic = substantive_num,
                              topic_name = substantive_top)

topics_and_nums$category <- ""
topics_and_nums$category[topics_and_nums$topic %in% c(9,28,39,63,65,71)] <- "Budgetary Politics"
topics_and_nums$category[topics_and_nums$topic %in% c(4,13,54)] <- "Civil Rights"
topics_and_nums$category[topics_and_nums$topic %in% c(1,11,15,16,42)] <- "Environmental Politics"
topics_and_nums$category[topics_and_nums$topic %in% c(17,26,27,34,48,58,73)] <- "IR"
topics_and_nums$category[topics_and_nums$topic %in% c(19,59,61)] <- "Political Institutions"
topics_and_nums$category[topics_and_nums$topic %in% c(23,29,44,49,62)] <- "Regulatory Politics"
topics_and_nums$category[topics_and_nums$topic %in% c(3,18,24,46,53,55,56,68)] <- "Welfare"
topics_and_nums$proportions <- colMeans(stm74$theta[,topics_and_nums$topic])

plot.STM(stm74,n = 5
         ,topics = topics_and_nums$topic[topics_and_nums$category=="Welfare"]
         , type = "summary", labeltype = "frex"
         ,topic.names = paste(topics_and_nums$topic_name[topics_and_nums$category=="Welfare"],"",sep = ":")
         ,main = "Welfare"
         ,xlim = c(0,0.08)
)

#####INDIVIDUAL CHARACTERISTICS#####
#####PARTISANSHIP#####
partyEffects <- extract.estimateEffect(prep, covariate = "Party", model = stm74,
                                       method = "difference", 
                                       cov.value1 = 1, cov.value2 = 0,ci.level = (1-0.05/74)
)
partyEffects <- filter(partyEffects,partyEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
partyEffects$label <- as.character(partyEffects$label)
partyEffects <- left_join(partyEffects, topics_and_nums, by = c("topic"))


plot(prep, covariate = "Party", 
     topics = topics_and_nums$topic[topics_and_nums$category=="Welfare"],
     model = stm74, method = "difference",
     cov.value1 = 1, cov.value2 = 0,
     ci.level = (1-0.05/74),
     xlab = "<-- Democrats --------------------------------------------- Republicans -->",
     main = "Effect of Partisanship on Topic Choice",
     xlim = c(-.01, .01),
     labeltype = "custom",
     custom.labels = topics_and_nums$topic_name[topics_and_nums$category=="Welfare"],
     verbose.labels = F, cex = 1)


#Plotting
partyEffects %>%
  ggplot(
    aes(x = estimate , y = reorder(topic_name,topic,decreasing = T)
        ,xmin = ci.lower, xmax = ci.upper
    ))+
  geom_point(fill = "black")+
  geom_errorbarh(aes(xmax=ci.upper, xmin=ci.lower),colour='black')+
  labs(x = "<-- Democrats ........................................................................................................ Republicans -->"
       ,y = "Expected Topic Proportion"
       ,title = "Effect of Party Affiliation on Topic Prevalence"
  ) +
  xlim(-0.01,0.01)+
  geom_vline(xintercept = 0, col='black',size=0.75)+
  theme(legend.position = "bottom")

#####SEX#####
sexEffects <- extract.estimateEffect(prep, covariate = "gender", model = stm74,
                                       method = "difference", 
                                       cov.value1 = 1, cov.value2 = 0, ci.level = (1-0.05/74)
)
sexEffects <- filter(sexEffects,sexEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
sexEffects$label <- as.character(sexEffects$label)
sexEffects <- left_join(sexEffects, topics_and_nums, by = c("topic"))


plot(prep, covariate = "gender", 
     topics = topics_and_nums$topic[topics_and_nums$category=="Welfare"],
     model = stm74, method = "difference",
     cov.value1 = 1, cov.value2 = 0,
     ci.level = (1-0.05/74),
     xlab = "<-- Female --------------------------------------------- Male -->",
     main = "Effect of Sex on Topic Choice",
     xlim = c(-.02, .02),
     labeltype = "custom",
     custom.labels = topics_and_nums$topic_name[topics_and_nums$category=="Welfare"],
     verbose.labels = F, cex = 1)

sexEffects %>%
  ggplot(
    aes(x = estimate , y = reorder(topic_name,topic,decreasing = T)
        ,xmin = ci.lower, xmax = ci.upper
    ))+
  geom_point(fill = "black")+
  geom_errorbarh(aes(xmax=ci.upper, xmin=ci.lower),colour='black')+
  labs(x = "<-- Female MCs ........................................................................................................ Male MCs -->"
       ,y = "Expected Topic Proportion"
       ,title = "Effect of Gender on Topic Prevalence"
  ) +
  xlim(-0.01,0.01)+
  geom_vline(xintercept = 0, col='black',size=0.75)+
  theme(legend.position = "bottom")

#####Race#####
raceEffects <- extract.estimateEffect(prep, covariate = "black", model = stm74,
                                     method = "difference", 
                                     cov.value1 = 1, cov.value2 = 0
)
raceEffects <- filter(raceEffects,raceEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
raceEffects$label <- as.character(raceEffects$label)
raceEffects <- left_join(raceEffects, topics_and_nums, by = c("topic"))


plot(prep, covariate = "black", 
     topics = topics_and_nums$topic[topics_and_nums$category=="Welfare"],
     model = stm74, method = "difference",
     cov.value1 = 1, cov.value2 = 0,
     ci.level = (1-0.05/74),
     xlab = "<-- NonBlack --------------------------------------------- Black -->",
     main = "Effect of Race on Topic Choice",
     xlim = c(-.02, .02),
     labeltype = "custom",
     custom.labels = topics_and_nums$topic_name[topics_and_nums$category=="Welfare"],
     verbose.labels = F, cex = 1)

#####Hispanic#####
HispanicEffects <- extract.estimateEffect(prep, covariate = "hispanic", model = stm74,
                                     method = "difference", 
                                     cov.value1 = 1, cov.value2 = 0
)
HispanicEffects <- filter(HispanicEffects,HispanicEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
HispanicEffects$label <- as.character(HispanicEffects$label)
HispanicEffects <- left_join(HispanicEffects, topics_and_nums, by = c("topic"))


plot(prep, covariate = "hispanic", 
     topics = topics_and_nums$topic[topics_and_nums$category=="Welfare"],
     model = stm74, method = "difference",
     cov.value1 = 1, cov.value2 = 0,
     ci.level = 0.9993243,
     xlab = "<-- NonHispanic --------------------------------------------- Hispanic -->",
     main = "Effect of Hispanic Ethnicity on Topic Choice",
     xlim = c(-.02, .02),
     labeltype = "custom",
     custom.labels = topics_and_nums$topic_name[topics_and_nums$category=="Welfare"],
     verbose.labels = F, cex = 1)

#####AGE#####
ageEffects <- extract.estimateEffect(prep, covariate = "age", 
                                      model = stm74,
                                     ci.level = 1-(0.05/74),
                                      method = "continuous")
ageEffects <- filter(ageEffects,
                      ageEffects$topic %in% 
                        topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ageEffects$label <- as.character(ageEffects$label)
ageEffects <- left_join(ageEffects, topics_and_nums, by = c("topic"))

plot(prep, covariate = "age", topics = topics_and_nums$topic[
  topics_and_nums$category=="Welfare"],
     model = stm74, method = "continuous",
     ci.level = 1-(0.05/74),
     npoints = 1000,
     xlab = "MC Age",
     main = "Effects of Age",
     #ylim = c(0, .035),
     verbose.labels = F,
     labeltype = "custom",
     custom.labels = topics_and_nums$topic_name[
       topics_and_nums$category=="Welfare"])


ggplot(ageEffects, 
       aes(x = covariate.value, y = estimate,
       ymin = ci.lower, ymax = ci.upper,
       group=factor(topic_name),
       fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line()+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "MC Age",
       y = "Expected Topic Proportion",
       title="Welfare",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")


#####INSTITUTIONAL-CONTEXTUAL#####
#####Change over Time#####
start.date <- as.Date(min(meta$date),'%Y-%m-%d')
end.date <- as.Date(max(meta$date),'%Y-%m-%d')
date_seq100 <- seq.Date(start.date,end.date, length.out=100)

#Correct days variable
prep$data$days <- as.numeric(prep$data$days)
#Extract days effects
daysEffects <- extract.estimateEffect(prep, covariate = "days", 
                                      model = stm74,
                                      method = "continuous",ci.level = (1-0.05/74))
#Filter based on relevant topics
daysEffects <- filter(daysEffects,
        daysEffects$topic %in% 
          topics_and_nums$topic[topics_and_nums$category=="Welfare"])
daysEffects$label <- as.character(daysEffects$label)
daysEffects <- left_join(daysEffects, topics_and_nums, by = c("topic"))
daysEffects$date <- date_seq100

ggplot(daysEffects[daysEffects$category=="Welfare",], aes(x = date, y = estimate,
                                                          ymin = ci.lower, ymax = ci.upper,
                                                          fill=factor(covariate)))+
  facet_wrap(~ reorder(topic_name,topic))+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .5, show.legend = F)+
  geom_line()+
  labs(x = "Time (1994-2021)",
       y = "Expected Topic Proportion"
  ) +
  labs(x = "Time (1994-2021)",
       y = "Expected Topic Proportion",
       title="Effect of Time on Topic Prevalence"
  ) +
  scale_x_date(date_breaks = "2 years",date_labels = '%y')+
  theme(legend.position = "bottom")

#####PARTY CONTROL#####
partyControlEffects <- extract.estimateEffect(prep, covariate = "partyControl", 
                                              model = stm74,method = "difference",
                                              cov.value1 = 1, cov.value2 = 0,ci.level = 1-(0.05/74))
partyControlEffects <- filter(partyControlEffects,partyControlEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
partyControlEffects$label <- as.character(partyControlEffects$label)
partyControlEffects <- left_join(partyControlEffects, topics_and_nums, by = c("topic"))


ggplot(partyControlEffects 
       ,aes(x = estimate , y = topic_name 
            ,xmin = ci.lower, xmax = ci.upper
       ))+
  geom_point(fill = "black")+
  geom_errorbarh(aes(xmax=ci.upper, xmin=ci.lower),colour='black')+
  labs(x = "Democratic Control ... Republican Control"
       ,y = "Expected Topic Proportion"
       ,title = "Effect of Party Control"
  ) +
  geom_vline(xintercept = 0, col='black')+
  theme(legend.position = "bottom")

plot(prep, covariate = "partyControl", 
     topics = topics_and_nums$topic[topics_and_nums$category=="Welfare"],
     model = stm74, method = "difference",
     cov.value1 = 1, cov.value2 = 0,
     xlab = "Democratic Control ... Republican Control",
     main = "Effect of Party Control",
     xlim = c(-.03, .03),
     labeltype = "custom",custom.labels = topics_and_nums$topic_name[topics_and_nums$category=="Welfare"],
     verbose.labels = F, cex = 1, ci.level = 1-(0.05/74))

#####Effects of Tenure#####
tenureEffects <- extract.estimateEffect(prep, covariate = "numberTerms", model = stm74,
                                        method = "continuous",
                                        ci.level = 1-(0.05/74))
#filter
tenureEffects <- filter(tenureEffects,
                     tenureEffects$topic %in% 
                       topics_and_nums$topic[topics_and_nums$category=="Welfare"])
tenureEffects$label <- as.character(tenureEffects$label)
tenureEffects <- left_join(tenureEffects, topics_and_nums, by = c("topic"))

#plot
tenureEffects %>%
  ggplot(aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ reorder(topic_name,topic))+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line()+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "# of Terms",
       y = "Expected Topic Proportion",
       title="Effect of Tenure on Topic Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####Effects of Region#####
regionEffects <- extract.estimateEffect(prep, covariate = "region", model = stm74,
                                        method = "pointestimate",
                                        ci.level = 1-(0.05/74))
#filter
regionEffects <- filter(regionEffects,
                        regionEffects$topic %in% 
                          topics_and_nums$topic[topics_and_nums$category=="Welfare"])
regionEffects$label <- as.character(regionEffects$label)
regionEffects <- left_join(regionEffects, topics_and_nums, by = c("topic"))

#plot
plot(prep, covariate = "region", 
     topics = topics_and_nums$topic[topics_and_nums$category=="Welfare"][1],
     model = stm74, method = "pointestimate",
     xlab = "Region",
     main = "Effect of Geography",
     xlim = c(-.03, .03),
     verbose.labels = T, cex = 1, ci.level = 1-(0.05/74))


####
ggplot(regionEffects 
      ,aes(y= covariate.value,x=estimate 
           ,xmin = ci.lower, xmax = ci.upper,
           group=covariate.value,
           colour=covariate.value
      ))+
  geom_point(fill = "black")+
  geom_errorbarh(aes(xmax=ci.upper,
                     xmin=ci.lower))+
  facet_wrap(~topic_name,ncol=2)+
  labs(y = element_blank()
       ,x = "Expected Topic Proportion"
       ,title = "Effect of Geography"
       ,colour = "Region"
  ) +
  geom_vline(xintercept = 0, col='black')+
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        legend.position = "bottom")

#####CONSTITUENCY#####
#####%WHITE#####
ConstitWhiteEffects <- extract.estimateEffect(prep, 
      covariate = "prcntWhite", model = stm74,
      method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitWhiteEffects <- filter(ConstitWhiteEffects, ConstitWhiteEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitWhiteEffects$label <- as.character(ConstitWhiteEffects$label)
ConstitWhiteEffects <- left_join(ConstitWhiteEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitWhiteEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0.5)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % White",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Whiteness on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")


#####%Black#####
ConstitBlackEffects <- extract.estimateEffect(prep, 
                                              covariate = "prcntBlack", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitBlackEffects <- filter(ConstitBlackEffects, ConstitBlackEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitBlackEffects$label <- as.character(ConstitBlackEffects$label)
ConstitBlackEffects <- left_join(ConstitBlackEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitBlackEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0.5)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % Black",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Blackness on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####%Asian#####
ConstitAsianEffects <- extract.estimateEffect(prep, 
                                              covariate = "prcntAsian", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitAsianEffects <- filter(ConstitAsianEffects, ConstitAsianEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitAsianEffects$label <- as.character(ConstitAsianEffects$label)
ConstitAsianEffects <- left_join(ConstitAsianEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitAsianEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0.5)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % Asian",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Asianness on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####%Hispanic#####
ConstitHispanicEffects <- extract.estimateEffect(prep, 
                                              covariate = "prcntHisp", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitHispanicEffects <- filter(ConstitHispanicEffects, ConstitHispanicEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitHispanicEffects$label <- as.character(ConstitHispanicEffects$label)
ConstitHispanicEffects <- left_join(ConstitHispanicEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitHispanicEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0.5)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % Hispanic",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Hispanicness on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####%BA#####
ConstitBAEffects <- extract.estimateEffect(prep, 
                                              covariate = "prcntBA", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitBAEffects <- filter(ConstitBAEffects, ConstitBAEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitBAEffects$label <- as.character(ConstitBAEffects$label)
ConstitBAEffects <- left_join(ConstitBAEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitBAEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0.5)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % BA",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Education on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####%HS#####
ConstitHSEffects <- extract.estimateEffect(prep, 
                                              covariate = "prcntHS", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitHSEffects <- filter(ConstitHSEffects, ConstitHSEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitHSEffects$label <- as.character(ConstitHSEffects$label)
ConstitHSEffects <- left_join(ConstitHSEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitHSEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0.5)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % HS",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Education on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####%Old#####
ConstitOldEffects <- extract.estimateEffect(prep, 
                                              covariate = "prcntOld", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitOldEffects <- filter(ConstitOldEffects, ConstitOldEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitOldEffects$label <- as.character(ConstitOldEffects$label)
ConstitOldEffects <- left_join(ConstitOldEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitOldEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0.5)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % Old",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Elderly on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####%Unemp#####
ConstitUnempEffects <- extract.estimateEffect(prep, 
                                              covariate = "prcntUnemp", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitUnempEffects <- filter(ConstitUnempEffects, ConstitUnempEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitUnempEffects$label <- as.character(ConstitUnempEffects$label)
ConstitUnempEffects <- left_join(ConstitUnempEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitUnempEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % Unemp",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Unemployment on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####%Recent#####
ConstitRecentEffects <- extract.estimateEffect(prep, 
                                              covariate = "recentArrivalPrcnt", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitRecentEffects <- filter(ConstitRecentEffects, ConstitRecentEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitRecentEffects$label <- as.character(ConstitRecentEffects$label)
ConstitRecentEffects <- left_join(ConstitRecentEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitRecentEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency % Recent Arrival",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Turnover on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")

#####totalEmploymentPop#####
ConstitTotalEmploymentPopEffects <- extract.estimateEffect(prep, 
                                              covariate = "totalEmploymentPop", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitTotalEmploymentPopEffects <- filter(ConstitTotalEmploymentPopEffects, ConstitTotalEmploymentPopEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitTotalEmploymentPopEffects$label <- as.character(ConstitTotalEmploymentPopEffects$label)
ConstitTotalEmploymentPopEffects <- left_join(ConstitTotalEmploymentPopEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitTotalEmploymentPopEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency Total Employment Pop.",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Population on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")



#####%medianIncome#####
ConstitmedianIncomeEffects <- extract.estimateEffect(prep, 
                                              covariate = "medianIncome", model = stm74,
                                              method = "continuous", nsims = 100,ci.level = 1-(0.05/74)
)
ConstitmedianIncomeEffects <- filter(ConstitmedianIncomeEffects, ConstitmedianIncomeEffects$topic %in% topics_and_nums$topic[topics_and_nums$category=="Welfare"])
ConstitmedianIncomeEffects$label <- as.character(ConstitmedianIncomeEffects$label)
ConstitmedianIncomeEffects <- left_join(ConstitmedianIncomeEffects, topics_and_nums, by = c("topic"))

#Plot
ggplot(ConstitmedianIncomeEffects, 
       aes(x = covariate.value, y = estimate,
           ymin = ci.lower, ymax = ci.upper,
           group=factor(topic_name),
           fill=factor(topic_name)
       ))+
  facet_wrap(~ topic_name)+
  #facet_grid(test~topic_name,margins = T,)+
  geom_ribbon(alpha = .4, show.legend = F)+
  geom_line(aes(color=factor(topic_name)))+
  geom_hline(yintercept = 0)+
  #geom_smooth(aes(colour=factor(topic_name)))+
  labs(x = "Constituency Median Income",
       y = "Expected Topic Proportion",
       title="Effect of Constituency Income on Welfare Prevalence",
       fill="Topic Name",
       colour="Topic Name"
  ) +
  theme(legend.position = "bottom")


#####Regression Tables#####
prep_summary <- summary(prep,topics=topics_and_nums$topic[topics_and_nums$category=="Welfare"])

print.summary.estimateEffectNames <- function(x, topic_names, digits = max(3L, getOption("digits") - 3L), 
                                              signif.stars = getOption("show.signif.stars"), ...) {
  cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"), 
      "\n\n", sep = "")
  
  for(i in 1:length(x$tables)) {
    cat(sprintf("\nTopic %i: %s", x$topics[i],
                topic_names[i]))
    cat("\nCoefficients:\n")
    coefs <- x$tables[[i]]
    stats::printCoefmat(coefs, digits = digits, signif.stars = signif.stars, 
                        na.print = "NA", ...)
    cat("\n")
  }
  invisible(x)
}

print.summary.estimateEffectNames(prep_summary,
                                  topic_names = topics_and_nums$topic_name[topics_and_nums$category=="Welfare"])

prep_summary_all <- summary(prep,topics=topics_and_nums$topic[1:18])
print.summary.estimateEffectNames(prep_summary_all, topics = topics_and_nums$topic[1:18],
                                  topic_names = topics_and_nums$topic_name[1:18])

prep_summary_all <- summary(prep,topics=topics_and_nums$topic[19:37])
print.summary.estimateEffectNames(prep_summary_all, topics = topics_and_nums$topic[19:37],
                                  topic_names = topics_and_nums$topic_name[19:37])